Towards a large-scale quantum simulator on diamond surface at room temperature 
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Strongly-correlated quantum many-body systems exhibits a variety of exotic phases with long- 
range quantum correlations, such as spin liquids and supersolids. Despite the rapid increase in 
computational power of modern computers, the numerical simulation of these complex systems be- 
comes intractable even for a few dozens of particles. Feynman's idea of quantum simulators offers 
an innovative way to bypass this computational barrier. However, the proposed realizations of such 
devices either require very low temperatures (ultracold gases in optical lattices, trapped ions, super- 
conducting devices) and considerable technological effort, or are extremely hard to scale in practice 
(NMR, linear optics). In this work, we propose a new architecture for a scalable quantum simulator 
that can operate at room temperature. It consists of strongly-interacting nuclear spins attached 
to the diamond surface by its direct chemical treatment, or by means of a functionalized graphene 
sheet. The initialization, control and read-out of this quantum simulator can be accomplished with 
nitrogen- vacancy centers implanted in diamond. The system can be engineered to simulate a wide 
variety of interesting strongly-correlated models with long-range dipole-dipole interactions. Due to 
the superior coherence time of nuclear spins and nitrogen- vacancy centers in diamond, our proposal 
offers new opportunities towards large-scale quantum simulation at room temperatures. 



Many intriguing phenomena in condensed-matter sys- 
tems originate from the interplay of strong interactions 
and frustrations. A representative example is frustrated 
quantum magnetism, where the spins cannot order to 
minimize all local interactions, and the ground state is 
highly degenerate [1, 2]. Together with long-range in- 
teractions between non-nearest neighbors, the frustrated 
quantum models give rise to intriguing quantum phases, 
e.g. supersolid [3, 4]. Moreover, they can also stabilize 
the long- sought quantum spin liquid [5-7] , which has con- 
nections with high-temperature superconductivity [8]. 
However, the properties of these systems have proven to 
be very hard to understand from numerical calculations, 
partly due to the combination of long-range quantum cor- 
relations and the superposition principle of quantum me- 
chanics [9] . This principle implies that the required com- 
putational resources grow exponentially with the num- 
ber of particles, making numerical approaches inefficient. 
Richard Feynman's idea [10] of quantum simulation pro- 
vides a powerful solution to this problem: one could gain 
deep insight into complex states of matter by experimen- 
tally simulating them with another well-controlled quan- 
tum system [11]. Large scale quantum simulations is ex- 
pected to be a powerful tool [12] for the investigation of 
fundamental problems in condensed-matter physics. 

Quantum simulation has attracted extensive research 
interest in the last decade. Various architectures for 
quantum simulations have been constructed based on dif- 
ferent systems, ranging from ultracold neutral atoms [13- 
15], trapped ions [16, 17], and photonic systems [18] to 
superconducting circuits in solid-state devices [19]. The 
still challenging goal is to realize a large-scale quantum 
simulator that cannot be efficiently solved numerically 



with classical computers, which will probably require go- 
ing beyond one-dimensional systems. In this work, we 
propose a scalable architecture that is of practical in- 
terest for large-scale quantum simulation. Our quan- 
tum simulator is based on lattices of interacting nu- 
clear spins, which can be fabricated chemically on the 
diamond surface [20, 21] or by depositing functional- 
ized graphene films [22-24]. We propose to implant the 
nitrogen- vacancy (NV) centers [25] beneath the diamond 
surface, and exploit them as an efficient control element 
for cooling, spin-spin interaction engineering, and read- 
out of the nuclear-spin quantum simulator. 

We explain how this simulator is constructed, and es- 
tablish schemes for its initialization, control and read- 
out. We analyze its validity by detailed numerical stud- 
ies, thus demonstrating the feasibility of our proposal 
within the current experimental capabilities. 

Construction of Hardware- We discuss two main paths 
for the fabrication of the hardware for our quantum sim- 
ulator. Firstly, large-scale lattices of nuclear spins can 
be constructed by chemically-controlled termination of 
diamond surfaces. The fluorine (^^F with spin ^), oxy- 
gen (^^O with spin |) and hydrogen/ hydroxy 1 group (^H 
with spin ^ and ^H with spin 1) termination of the dia- 
mond surface can be obtained from the process of chem- 
ical vapor deposition (CVD) ,or by functionalizing the 
diamond surface. As a representative example, we will 
mainly concentrate on fluorine-terminated diamond sur- 
face, which has a positive electron affinity. The two most 
important diamond surfaces are the (111) and (100) sur- 
faces, which constitute the crystal faces of polycrystalline 
CVD diamond films and can be selectively grown with 
appropriately controlled process parameters [20]. The 



FIG. 1: Lattices of fluorine nuclear spin, (a) A triangular nuclear spin lattice on the fluorine-terminated diamond (111) 
surface, (b) A rectangular nuclear spin lattice on the fluorine-terminated diamond (100) 2x1 surface. The distance between 
two nearest neighbour fluorine atoms is about 2.5 (unit A), (c) A double-layer triangular lattice from fluorographene. Yellow 
atoms represent carbon, and green atoms represents fluorine. 



(Ill) surface of diamond is the natural cleaving plane 
of diamond, and has one dangling bond per surface car- 
bon atom which is terminated by carbon-fluorine bonds 
on the fluorine-terminated diamond surface [26]. The 
fluorine atoms naturally arrange in a triangular lattice 
with nearest-neighbour distances of about 2.5A [20, 21], 
see Fig. 1(a). The (100) surface of diamond shows two 
dangling bonds per surface carbon atom, which will un- 
dergo a reconstruction into a 2 x 1 geometry with neigh- 
boring surface carbon atoms forming a 7r-bonded dimer. 
The remaining dangling bonds are terminated by carbon- 
fluorine bonds, which leads to a rectangular lattice of 
fluorine atoms [20, 21], see Fig. 1(b). Functionalized 
graphene (fluorographene) provides a double-layer trian- 
gular lattice of fluorine atoms, see Fig. 1(c). This can be 
obtained through the mechanical cleavage of graphite flu- 
oride, or by exposing graphene to atomic fluorine formed 
by decomposition of xenon difluoride (XeF2) [23]. 

In addition to the nuclear spins, we shall introduce NV 
centers by shallow implantation a few nanometers below 
the surface of diamond [27, 28]. This constitutes a fun- 
damental ingredient of our quantum simulator, since it 
allows for the initialization and read-out of the nuclear 
spins. Let us remark that in contrast to graphene, fluoro- 
graphene exhibits a large band gap of more than 3.5eV, 
which is larger than the optical gap of NV centers in dia- 
mond (^ 2.9eV). This avoids the unwanted fluorescence 
quenching of NV centers and is thus crucial for using the 
NV centers for the control of the nuclear spin arrays. 

Engineering of Interacting Hamiltonian- The nuclear 
spins interact with each other via magnetic dipole-dipole 
interactions as Vij = g{rij)[si ■ Sj — 3 (s^ • Vij) {sj ■ rij)], 
where are the nuclear spin operators, g{rij) = 
(fi^/io7i7i)/(47rrf^-), 7^, 7^ are gyromagnetic ratios of the 
ith and jth nuclear spin, which are connected by the 
vector Vij = rijVij. Since nearest-neighbor fluorine nu- 
clear spins are at a distance of 2. 5 A, a coupling strength 
of g{a) 6.8kHz is achieved. A strong magnetic field, 
which leads to energy level shifts exceeding the nuclear 
spin coupling strength, simplifies the effective Hamilto- 
nian due to the rotating wave approximation (RWA) to 



an XXZ model 

y,, = g{r,,) [s^s| - A (s^s^ + s^)] (1) 

with A = |. Our calculations verify that such a RWA 
can already be well satisfied for a magnetic field as low as 
750G (3MHz), see SI for details. We denote the diamond 
surface, on which the nuclear spin lattice is constructed 
as the X-Y plane, while the vector perpendicular to it 
defines the Z axis (spatial axes), and the magnetic field 
direction as rh axis which gives the quantization of nu- 
clear spins. The spin-spin interaction strength can be 
tuned by changing the direction of magnetic field as 

di^ij) = di^j) (1 - 3cos^%) (2) 

where 6ij is the relative angle between Vij and rh. Thus, 
by changing the direction of the external magnetic fields, 
we can control the spatial anisotropy and the sign of the 
interaction strength, which determines whether the inter- 
action is ferromagnetic or anti-ferromagnetic and thereby 
induces spin frustration. The value of the spin anisotropy 
A can be tuned with gradient magnetic fields to modu- 
late the hopping coupling (sfs^ + sj^sj) while keeping the 
repulsive interaction sfs| unchanged [29]. In this way the 
interplay between geometrical frustration and the effects 
of quantum fluctuations on the realisation of spin liquid 
phase could be tested [30] . As an alternative scheme with 
higher nuclear spin species, for example ^H with spin 1 




FIG. 2: Tuning spin anisotropy with dressed states. 

Simulation of an effective spin-| by the dressed states from 
a higher spin with radio frequency driving fields (with Rabi 
frequency Qh and detuning Ah). 
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FIG. 3: Initialization and read-out of nuclear spins. 

NV center lying beneath the nuclear spin layer in bulk dia- 
mond is exploited to initialize and measure nuclear spins on 
diamond surface. The spin-1 ground state of NV center can 
be optically readout via spin-dependent fluorescence, and can 
be coherently controlled with microwave fields. 

and ^''O with spin-|, one can tune the spin anisotropy 
A by applying continuous fields corresponding to the nu- 
clear spin transitions |m — 1) ^ |m) and |m) ^ |m + 1) 
with Rabi frequencies Qh and detuning A/^, see Fig. 2. 
An effective spin- ^ can be encoded in two of the induced 
dressed states with the spin anisotropy A flexibly tuned 
(see SI for details). We remark that spin Hamiltonians 
for higher spins may directly lead to rich quantum phases, 
see e.g. [31, 32]. 

Cooling of Quantum Simulator- The reliable prepa- 
ration of the quantum simulator in a low entropy state 
is a prerequisite for the observation of quantum phases. 
Here we explain in detail how the nuclear spin lattice can 
be initialized to a well defined spin direction using dy- 
namical nuclear polarization via NV centers in diamond, 
see Fig. 3. The entire process consists of repetitive cy- 
cles. In each cycle, the NV center is first prepared in the 

1+) = y^(|ms = 0) + \ms = +1)) state by optical pump- 
ing to the \ms = 0) state using a green laser (532nm) and 
a subsequent | microwave pulse. A microwave field (with 
Rabi frequency ^nv) on resonance with the electronic 
transition \ms = 0) ^ \ms = +1) will lock the NV elec- 
tron spin state. If the Hartmann-Hahn condition with 
nuclear spins [33] is satisfied, the NV electron spin polar- 
ization will be transferred to the nuclear spins [34] (see 
SI for details). Due to their proximity, the interactions 
between nuclear spins exceeds that between the electron 
spin of the NV center and the nuclear spins. Therefore, it 
will be advantageous to effectively decouple nuclear spin 
interactions, which also facilitates estimates of the cool- 
ing efficiency. To this end, we apply a radio frequency 
(RF) field with the amplitude Qp whose frequency is de- 
tuned from the Larmor frequency of the nuclear spins 
by Ap. The eigenstates of the effective local Hamilto- 
nian of nuclear spin Hi = (l]p/2) sf + Apsf gives a new 
spin basis Written in such a spin basis, the 



energy conserving terms of nuclear spin interactions can- 
cel each other when ftp = 2V^Ap, and the anti-rotating 
terms are suppressed as long as the energy mismatch = 
(l]p + Ap)^ Qij is fulfilled. Thus, the nuclear spins be- 
have as isolated particles, and the Hartmann-Hahn con- 
dition becomes ^nv = (^nB - Ap) + [(l]p/2)^ + ^p]^, 
see SI for more details. Such a mechanism for the iso- 
lation of the nuclear spins (as we will discuss in the fol- 
lowing section) can also reduce the perturbation on the 
nuclear spin state during the read-out process, which will 
be beneficial for the accuracy of the measurement. 

To verify the validity of this idea, we have used a 
Chebyshev expansion [35] to calculate polarization dy- 
namics with the exact Hamiltonian of a 3 x 3 nuclear spin 
lattice assuming a distance of 5nm from the NV cen- 
ter. It can be seen from Fig. 4(a) that the coupling be- 
tween nuclear spins is effectively eliminated. We compare 
the exact numerical calculation with the results for iso- 
lated nuclear spins under the spin temperature approx- 
imation in which coherences between nuclear spins are 
neglected [36], see Fig. 4(b), which show a good agree- 
ment. To achieve an ultra-high polarization given by the 
spin temperature approximation, one can introduce mag- 
netic noise to remove coherence among nuclear spins in 
between polarization cycles (during which the NV cen- 
ter is polarized to the state \ms =0) and would not 
affect nuclear spins). From the spin temperature ap- 
proximation, one can estimate the polarization rate as 

Pi = {gifr = eigif/^T^kigiY « U^gi), where T is 
the polarization cycle time, e is a chosen constant, and 
is the flip-flop rate between the NV electron spin and 
the nuclear spin k (see SI for details). The required po- 
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FIG. 4: Isolating and initialization of nuclear spins, (a) 

The localization (circle, with RF-driving) and de-localization 
(square, without RF-driving) of spin excitation initially cre- 
ated in one specific site, demonstrating the efficiency of decou- 
pling of nuclear- nuclear interaction in the former case. The 
driving parameter is uj^ — 160kHz. (b) The average nuclear 
spin down state population as a function of the polariza- 
tion time t: Exact numerical simulation (purple, square) spin 
temperature approximation (yellow, diamond) with the cy- 
cle time T = 30/is; Exact numerical simulation (red, circle) 
with r = 120/is and magnetic noise is introduced to remove 
nuclear spin coherence every 20 cycles. The nuclear spin lat- 
tice is 3 X 3 with a distance of 5nm to the NV center. The 

magnetic field direction is rh — (\/|, 0, \ \)- 
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FIG. 5: Quantum magnetic phase transitions of the fluorine quantum simulation on a triangular lattice, (a) 

Tuning from ID anti- ferromagnetic chains to ferromagnetic order with the magnetic field in the direction m — (0, cos^, sin^). 
(b) The phase transition from ferromagnetic- antiferromagnetic (F-AF) order to ferromagnetic- (alternative) antiferromagnetic 
(F-NAF) controlled by the magnetic field in the direction m = (cos^, 0, sin^). (al-bl) The spin structure factors are calculated 
using the Lanczos algorithm for a 6x4 nuclear spin lattice under periodic boundary condition. (a2-b2) Comparison of the spin 
structure factors from exact diagonalization of a 4 x 4 nuclear spin lattice (lines) and the estimated values (which have been 
rescaled for visibility) via NV measurements (circles and squares) by numerical simulations. The interaction time (between the 
electron spin of NV center and the nuclear spins) for read-out is r = 60/is. 



larization time scales linearly with the total number of 
nuclear spins N and the inverse effective temperature. 
The ultimate polarization efficiency will be limited by 
the relaxation time Ti of nuclei, which can be as long as 
a few hours even at room temperature. The polarization 
cycle time sets the required coherence time of NV centers 
and nuclear spins to a few hundreds of /is, which is read- 
ily achievable with the current experimental techniques in 
diamond samples. The polarization efficiency can also be 
improved by optimizing the polarization cycle time and 
exploiting several NV centers. Once the nuclear spins 
have been initialized, by performing an adiabatic pas- 
sage, one can prepare the system in the ground state of 
the engineered interacting Hamiltonian (similar to other 
types of quantum simulator [37]). 

Measurement of Quantum Simulator- Before dis- 
cussing the static and dynamical properties of the pro- 
posed quantum simulator, let us describe the measure- 
ment schemes and demonstrate their viability by means 
of numerical simulations. The main goal is to measure 
observables that can provide information about nuclear 
spin state [38] . It is challenging to measure nuclear spins 
directly due to their small magnetic moments. However, 
NV centers implanted beneath the diamond surface will 
provide the solution as a measurement interface for nu- 
clear spin states [34]. Before the measurement, we apply 
a RF-pulse to map the nuclear spin basis from |t) and 
1^) to It) and 1^), in which the nuclear spins can be effec- 
tively decoupled from each other by a continuous driving 
field as described in the process of initialization. The 
NV center is prepared in the state (/i = ±), and is 
driven to match the Hartmann-Hahn condition between 
the electron spin of NV center and the nuclear spins. Af- 
ter the electron spin interacts with the nuclear spins for 
time r, we measure the population of state of the NV 



center, which is given by = Ei Ej(^i^^/)(sf sp 
to second-order in r, where /i, = ±. The above ob- 
servables include both local contributions of individual 
nuclear spins (for i = j) and two-point correlations of nu- 
clear pairs (for i 7^ j). We can extract information about 
the average nuclear spin magnetization by Am = P- — 
= 2r^ '^Ziidi')'^ i^z) ^ transverse correlation as 

AxY = P++P+ = 2T2EiEi(5.^5,^)(s?s| + sfsp. The 
correlation function along the other directions Ayz and 
Axz can be obtained by applying the Hadamard oper- 
ation Oh = it.)(ti + i;.)ai with It,) = ^(it) + It)), 

and the phase transformation 0/ = |t)(t| +^|t)UI on the 
nuclear spin state before the measurement. To measure 
observables such as structure factors, which are very im- 
portant for the study of quantum phase transition and 
for inferring entanglement properties [39], we can use a 
gradient field, with which the nuclear spin at the po- 
sition experiences a magnetic field bi = Vi ■ ihx^hy) 
and gains a position-depend phase ipi 
q = jNtp{bx,by) after a certain time tp 
ing the same measurement as before, we can obtain 
A^y = t2 EMat) COS iivi - r,) ■ q){sfs^ + s^'sj) 
and Ay^, A^^ in a similar way. We remark that it 
is possible to generate a gradient field as large as 15G 
(about 10 times larger than the coupling between fluo- 
rine nuclear spins) over the lattice constant (2. 5 A) by 
one single magnetic tip with the state-of-art technology 
[40]. The validity of our measurement scheme is numeri- 
cally tested in the context of witnessing quantum phase 
transitions as described in the following section. 

Frustrated Quantum Magnetism - Our system can sim- 
ulate quantum spin models with the tunable spin-spin 
interaction g{Tij): positive g{Tij) correspond to anti- 
ferromagnetic (AF) interaction, and negative g{Tij) is 



yJ 

Ti • q where 
By perform- 
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ferromagnetic (F) interaction. In the triangular lattice 
of the fluorine simulator, the nearest-neighbor nuclear 
spin interactions are denoted as Ja^, Jas, Ja^ with ai = 
(1,0,0) a2 = (i^,0), as = (-^^,0). The long- 
range interactions and spin frustration make it hard to 
perform numerical simulations using the quantum Monte 
Carlo method due to the subtle sign problem [41]. Here, 
we exactly diagonalize the system on a 6x4 lattice using 
the Lanczos algorithm under periodic boundary condi- 
tion to provide evidence for various phases of quantum 
magnetism. 

In the situation where Ja^ = Ji = ga is positive 
(AF) , and Ja^ = Jds = J2 = ^a(l - | cos^ (9) is 
negative for cos^ > | (the magnetic field direction 
is m = (0, cos^, sin^)), the triangle which consists of 
J2 — J2 — Ji is spin frustrated, see Fig. 5(a). For small val- 
ues of I J2/ Ji I, the system consists of ID (AF) chains with 
weak intra-chain interaction which induces ferromagnetic 
order in the sublattice of every two ID chains and is 
characterized by the (normalized) spin structure fac- 
tor %(^,0) = {l/N^){{T.^^'^'^'^l^{T.i^'^'^'^i)^) with 
q = (7r,0). As | J2/J1I increases, the competition between 
anti- ferromagnetic (Ji) and ferromagnetic (J2) interac- 
tions lead to the ferromagnetic phase above the point 
\J2/Ji\ = 1 corresponding to the spin structure factor 
5'g*(0,0), see Fig.5(ai). Note that the largest non-nearest 
neighbor interaction J3 is ferromagnetic and is essential 
for the emergence of the ferromagnetic phase which is 
absent for the short-range model, see SI for details. 

For the other situation where Ja^ = Jas = J2 = 
g{a){l — I cos^ ^) is always positive (AF), and Ja^ = 
Ji = g{a){l — 3cos^0) is negative (F) for cos^ > |, 
where is defined by the magnetic field direction is m = 
(cos ^,0, sin 0) . Considering only the nearest-neighbor in- 
teractions, the system is non-frustrated, and is expected 
to be ferromagnetic in the ai, while anti- ferromagnetic in 
the a2 direction (F-AF phase), which is identified by the 
spin structure factor 5'g*(0, 27r/v^). As the value of cos^ 
approaches to 1, the non- nearest neighbor interactions 
J31 = 1/3a/3, J32 = (1 — |cos^^)/3\/3 become compa- 
rable to J2. The competition between them leads to a 
magnetic phase identified with the spin structure factor 
5'<f(0,7r/v^). The spins are ferromagnetic in the ai di- 
rection, while they are anti-ferromagnetic between next- 
nearest-neighbor chains in the a2 direction, see Fig.5(bi). 

We also check the feasibility of using NV centers to 
identify different magnetic phases. We numerically cal- 
culate the dynamics during measurement and obtain the 
estimation of spin structure factors. Due to the limit 
of computational overhead, we consider an example of 
a 4x4 nuclear spin lattice with periodic boundary con- 
ditions. By applying a gradient field we can generate 
the relative phases between different nuclear spins corre- 
sponding to the spin structure factor, which can then be 
estimated by the observable S^{qx^qy) cx Axz{Qx^Qy) + 
^yz{Qx^ Qy) — ^xy{Qx^ Qy) via NV centers. We find that 



the estimation is in good agreement with the results from 
exact diagonalization, and that it can witness quantum 
phase transitions between different magnetic orders, see 
Fig.5(a2-b2) and SI for more details. 

Quantum Superfluid and Supersolid- The nuclear spin 
Hamiltonian can be mapped to the hard-core boson 
model by the Holstein-Primakoff transformation as 



^ij + ^^4) + M X] (3) 



with rij 



f + ^ {rii = 0, 1). Here, the chemical poten- 
tial is /i = (jnB) — ^jg{^ij)^ the repulsive interaction 



is Vij = gi^ij)^ and the hopping is ti 



rg{rij). Our 



proposed system can therefore simulate the hard-core bo- 
son model, which demonstrates interesting phases such as 
(long-range off-diagonal order) superfluid, and moreover 
a supersolid phase characterized by both long-range off- 
diagonal and diagonal order. Note that quantum sim- 
ulation of similar models with polar molecules in opti- 
cal lattices have been proposed [42, 43] and numerically 
studied in [44, 45] for long range repulsive dipole interac- 
tion (Vij). Our system possesses long range interactions 
(both Vij and t^j), and thus offers a platform to inves- 
tigate rich phases of hard-core bosons with potentially 
novel features. Indeed, these models with frustration and 
long-range interactions pose considerable challenges to 
the classical numerical simulations due to the sign prob- 
lem for 2D systems. 

With the magnetic field along the direction m = 
cos6>a3 + sin6>(0,0, 1) = (-^ cos6>, ^ cos6>, sin6>), the 
nearest-neighbor interactions are Vq,^ = Va^ = Vi = 
ga(\ - f cos^^), and V^, =¥2= ga(l - 3cos2^). By 
changing the value of the magnetic direction angle ^, we 
can gradually tune the geometric frustration as quanti- 
fied by the ratio V2/V1. For cos^ G [0, \J\]-, the values 
of V for all interactions have the same sign (including 
long-range interactions), and one can simulate such a 
model with the directed loop algorithm in the stochas- 




FIG. 6: Quantum phases of superfluid (SF) and super- 
solid (SS). (a) The average filling (n) as a function of chem- 
ical potential /i — /io (in the unit of ^a), where /xo is the chem- 
ical potential corresponding to half-filling, (b) The superfluid 
density ps and the normalized structure factor Sg*(7r, tt/v^)- 
The system size is 12 x 12 and t/V — 0.2. The magnetic 

field direction is ?fi = cos 6a2, + sin 6Z with cos — ^J\. The 
temperature used in our simulation is T / Qa — 0.1. 
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tic series expansion representation of the ALPS library 
[46]. By comparison with the short-range model, it can 
bee seen that the long-range interaction significantly en- 
hances the superfluidity, see SI for more details. By tun- 
ing the anisotropy value, we can observe quantum phase 
transitions between solid (S), supersolid (SS) and super- 
fluid (SF), see Fig. 6. One can also use a gradient fleld 
to selectively tune hopping interactions. For example, 
with a gradient fleld that decreases along the direction 
Ab = A5(^, ^), the hopping interaction will be sup- 
pressed except in the direction as = (— ^, ^). 

Outlook — We have introduced a novel platform for 
quantum simulations that incorporates the properties of 
diamond surfaces with the characteristic of NV centres 
in diamond. With numerically tractable examples, we 
have demonstrated that important quantum phases can 
be observed with the proposed quantum simulator. Due 
to the flexibility of the method in engineering Hamilto- 
nians, quantum models with novel features will be re- 
alizable, which would require more exhaustive classical 
simulation for understanding the richness of the quan- 
tum phases. In addition to the static properties of the 
quantum phase transitions, the present quantum simula- 
tor is also capable of studying quantum many-body dy- 
namics, such as quantum quenches and the generation 
of multi-particle entanglement, see SI for a simple exam- 
ple. The prerequisite technologies for the implementation 
of such a quantum simulator, such as the techniques for 
charge state manipulation of NV centers in (surface func- 
tionalized) diamond [47], coherent control of NV centers, 
have been developed already [48-50]. The present pro- 
posal will stimulate the interest of material scientists to 
fabricate other candidate systems, in particular for the 
simulation of quasi 3D systems (e.g. a number of fluo- 
rographene sheets). The relevant tools in this work may 
also be beneflcial for further research in coherent control 
on surfaces/mono-layer fllms. 
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Supplementary Information 

Hamiltonian of NV centers and nuclear spins. — The system under discussion comprises the electron spin of the NV 
center and the nuclear spin array on the diamond surface. The dynamics of these parts and their mutual interaction 
is described by the Hamiltonian 

H = Hil + Hp+H^y_p, (1) 

where the individual terms are discussed below. The lower energy manifold of the NV center, made up of the electron 
spin states, is described by the Hamiltonian 

nil = DSl - 7eB • S (2) 

where D = 2.87GHz is the zero field splitting, 7e is the gyromagnetic ratio of the electron spin, and S is the spin-1 
operator. Here, H^y describes three electron spin states, see Fig.l, where the degeneracy of the states m^ = +1 and 
m^ = — 1 is lifted by a magnetic field. Two of the eigenstates Im^ = —1) = |t) and Im^ = 0) = |^) can thus serve as 
an effective spin-^ with the following Hamiltonian 

rjgS ^nv /ox 

Hnv = -^^z (3) 
where = |t)(t| ~ is the Pauli operator. To achieve the goals of initialization, control and measurement of 



637 nm 

NV probe 

— '"^ ^ 




FIG. 1: The sketch of the energy levels of NV center. The spin-1 ground state of NV center (|ms = 0), \ms = ±1)) can be 
optically readout via spin-dependent fluorescence, and can be coherently controlled with microwave fields. We continuously 
drive the one of its electronic transition |ms = —1) ^ |nis = 0), which provides an effective spin— | serving as a quantum probe 
for target nuclear spins and in the mean time is decoupled from other noise. 

the quantum simulator, while at the same time decoupling the electron spin of the NV center from other spin species 
[SI], we continuously drive the electron spin of the NV center with a microwave field on resonance with the electronic 
transition [m^ = —1) = |t) ^ [m^ = 0) = |^). The effective Hamiltonian of the NV center then becomes 

Hil = ^cx. (4) 

where ^Inv is the Rabi frequency of microwave driving and ax = |t)(i| + li)(t|7 see Fig.l. The Hamiltonian of the 
nuclear spin quantum simulator on the diamond surface is given by 

% = E ^^B^^' + £ E ?l • - 3 • ^^J) i^J ■ ^^J)] (5) 

where is the gyromagnetic ratio of nuclear spins, and Vij = rijVij is the vector that connects two nuclear spins. The 
nuclear spin is quantized along the direction of the magnetic field. The interaction between the NV center electron 
spin and the nuclear spins is as follows 

Hm-F = £ E ^ [S • - 3 (S • r) (s, • f)] . (6) 

i ^ 

The Hamiltonian for other species of nuclear spins, e.g ^^C and ^^N of NV center, takes a similar form. For an 
isotopically engineered diamond, the number of ^^C can be reduced to 0.01% [S3], and only very few ^^C will 
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affect NV centers. The Hamiltonian for ^^N of NV center has a similar form with ^^C with the hyperfine couphng 
Aq = 2.1MHz and Aq = 2.3MHz [S4]. Note that the use of continuous wave decouphng techniques [S2] can reduce 
the impact of noise from impurities such as ^^C or ^^N further rendering their impact neghgible. 

Hamiltonian engineering for the quantum simulator. — The Hamiltonian for nuclear spins on the diamond surface 
lattice is presented in Eq.(5). We apply a radio-frequency field with the amplitude and the frequency uj = jn^—ujf 
where cuf denotes the detuning. Therefore, the Hamiltonian of Eq.(5) can be rewritten as 



Mo Tat 
47r ^ rf. 



[si • s^- - 3 (si • Yij) {sj ■ i-ij)] + 2nF cos [(7ArB - ujF)t] ^ 



The effective Hamiltonian in an interaction picture with respect to Hq = ^^{jn'B — ^F)sf is given by 



i 



47r ^ 



(7) 

(8) 
(9) 



where 



H. = 



7^ 



(10) 



with the spin anisotropy A = ^. Eq.(lO) is the quantum simulator Hamiltonian, for which we can tune the nuclear 
spin interaction strength and the spin anisotropy as described in the main text. If the amplitude of the RF field 
is much larger than the nuclear spin interaction, the ground state of the Hamiltonian will be \X) = • • • |^^), 
which can be prepared with NV centers (more details will be discussed in the following sections). By adiabatically 
decreasing the amplitude Qf^ the system will end in the ground state of the new Hamiltonian Hs. 

Besides using gradient fields, in combination with higher nuclear spin species (e.g. ^H with spin 1 and ^^O with 
spin- 1), it is also possible to tune the value of the spin anisotropy A for an effective spin-^ Hamiltonian by applying 
continuous fields with appropriate Rabi frequencies and detunings. The quadrupole together with Zeeman splitting 
allow us to apply microwave fields to selectively drive the nuclear spin transitions |m — 1) ^ |m) and |m) ^ |m + 1) 
with the detuning and —Ah respectively, see Fig. 2 of the main text. The Rabi frequencies are the same and 
denoted as fth- In the subspace of {|m — l),|m),|m+l)}, the Hamiltonian in the interaction picture is 



Ho 



Its three eigenstates are 



\d-i) = 


1 


m - 


hi) 


\do) = 


1 

^' 


;|m4 


-i> 


\d+i) = 


1 


m - 


hi) 



Ah 
^Ih 

rih Ah 



1 + y/l + 8(n;,/A;,)2 

|m-l)), 

1 - yi + 8(!^;,/A;,)2 



(11) 



|m) + |m- 1) 



Im) + Im- 1) 



(12) 
(13) 
(14) 



If we choose the dressed states {\do), or |<i-i), |<^+i)} as above to represent an effective spin-^, and express the 

nuclear spin interaction in this subspace, the value of the spin anisotropy J^/j" = A (with and j" denotes the 
strength of coupling sf ■ -\-s^ ■ and sf • s| respectively) will depend on the ratio between the Rabi frequency and 
the detuning 0.^/ A^. In Fig. 2, with the example of ^^O, we see that the value of J^/j" can be fully tuned in the 
range of (— oo, oo) by changing the value of Qh/ A^. 
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FIG. 2: The spin anisotropy /J^^ = —A as a function of the ratio between the radio frequency field amplitude and the 
detuning Qh/^h- (a-b) The effective spin-| is represented by two of the dressed states {|c?o), |<i+i)}; (c-d) The effective spin-| 
is represented by from the original nuclear spin states {| — + + of^^O. 

Validation of the rotating wave approximation. — We have adopted the rotating wave approximation (RWA) to 
obtain the XXZ spin Hamiltonian. More specifically, the Hamiltonian describing the magnetic dipole-dipole interaction 

Hp = J2 iNBst + ^ ^ ^ [s, . s, - 3 (s . r,,) (s, • f.,)] (15) 
is approximated by the effective Hamiltonian as 

Hs = Y: 7NBsf + g ^ 2| (1 - 3(r|,.)^) K • - A(sf • sj + sf • s^] (16) 

The rotating wave approximation is valid if the magnetic field is much larger than the spin-spin coupling, i.e. 
7ArB ^ (moTat) / i^'^'^ij)' We have verified the rotating wave approximation by comparing the dynamics from the 
Hamiltonian of Eq.(15) and Eq.(16). With the initial state |?/^) = |t) (8) ||) (8) • • • (8) ||), we plot in Fig. 3 the state 




t(ms) t(ms) t(ms) 

FIG. 3: Dynamics of spin state population Pj'^'* for the initial state — |t) 8) |i) • • • Complete Hamiltonian in 
Eq.(15) (curves) vs. effective Hamiltonian in Eq.(16) (circles). The nuclear spin lattice is 3 x 3, the magnetic field direction 
is m = (0,0, 1) and its strength is 750G (corresponding to 3MHz for fluorine nuclear spin). Note that the scales in different 
panels have different order of magnitude. 
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population for each nuclear spin. One can see that the dynamics of the complete Hamiltonian agrees well with the 
one from the effective Hamiltonian. 

Isolating nuclear spins for cooling and measurement. — Due to the very small energy splitting of nuclear spins 
subjected to a magnetic field, it is very challenging to cool nuclear spins to a sufficiently low entropy state, such that 
they exhibit quantum properties. To overcome this problem, we propose to use the electron spin of NV centers to 
transfer polarization to the nuclear spins. Hence, we can take advantage of the fact that the NV center electron spin 
can be easily polarized even at room temperatures. Due to the large zero field splitting of the NV center electron 
spin, direct polarization transfer to nuclear spins is not possible due to the large energy mismatch. One possible way 
to bypass this problem is to continuously drive the NV electron spin to induce microwave dressed states {Itx)^ \ix)} 
(i.e. the eigenstates of the effective Hamiltonian in Eq.(4) realizing an effective spin— ^), the polarization of which 
can now efficiently be transferred to the nuclear spins when the Hartmann-Hahn condition is fulfilled (i.e. when the 
driving Rabi frequency matches the Larmor frequency of the nuclei) [SI, S7]. It is important to note that the coupling 
strength between nearest-neighbour fluorine nuclear spins is 6.8kHz, while the interaction between the NV center 
and the fluorine nucleus is much smaller (e.g. ~ 0.5kHz for a distance of 5nm). Therefore, it is rather complicated 
to achieve the Hartman-Hahn condition with each energy gap as required for the polarization exchange. The most 
promising solution consists of decoupling the nuclear spins from each other. This is also important for the read-out 
of quantum simulation result in order to avoid significant changes of nuclear spin states during the measurement. To 
achieve this goal, we apply a RF-field with the amplitude denoted as l^p whose frequency is off-resonant with the 
Larmor frequency of the nuclear spins by the detuning Ap. The total Hamiltonian is written as follows 

Hp =^a. (H^l) (17) 

2 

+ J2 7ivBsf + 1^ ^ 2f [s, . s,- - 3 (si ■ Vij) {sj ■ hj)] + cos [{j^B - Ap)t] ^ sf (Hp) (18) 

i i,j ij i 

+ ^(I + a,)(^^(aX+«ysf +«.sf) (-H"nv-F)- (19) 

i 

We choose the magnetic field such that the energy splitting of the nuclear spin states is much larger than the nuclear 
coupling strength. This allows us to simplify the above Hamiltonian and obtain 

Hf =^<x. {Uf,) (20) 



+ ^ E«-«^' + JH) ® E [«)' + («r)']'^'sr + ^-^^ (^nv-F) (22) 

i i 

where cjnv = ^nv - (7a/"B - Ap), = |tx)Ux| and cr^~^ = ||a^)(t x\ are the raising and lowering operator for the 
dressed spin-^ of the NV center, and the term \ ^Obz^\ gives the additional field on nuclear spins as created by the 

NV center. We note that the local part of nuclear spin Hamiltonian i^p = ^ • + Ap ^ • sf introduces a new 
nuclear spin basis which depends on the ratio Ap/l]p. In particular, if 

l^p = 2V^Ap (23) 



the new nuclear spin basis is 



with 



It) = cos(^)|t)+sin(y))|i) (24) 
\i) = sin(^)|t)-cos((p)|i) (25) 



cos(2y)) = yi and sin(2(^) = ^ ^. (26) 
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Written in such a basis, the excitation number conserving terms from sf sj and sf + cancel each other due 
to their opposite signs, see Eq.(21), additionahy the effects of the other non-energy conserving terms are suppressed 

by the large energy mismatch as long as uj^ = -^Ap + {Qp/2y is much larger the nuclear spin interaction. For the 

fluorine nuclear spins, RF-fields with an amplitude as strong as 200kHz are available, which are much larger than the 
nearest-neighbor nuclear spin interaction (i.e. 6.8kHz). The effective Hamiltonian therefore becomes 

<^ = + -f E + E ^1 + E {<) + h.c) (27) 

i i i 

where sf = |t)(t| ~ ^+ = |t)(i|7 s~ = |i)(t| the Pauli operator and raising (lowering) operator for 

nuclear spins in the effective spin basis. The coefficients ^| = | cos{2(p)af represent the additional field induced by 

1/2 

the NV center on the nuclear spins, g:^ = | [1 + cos{2(p)] [(af )^ + gives the rate of polarization exchange 

between the NV center and the nuclear spins, where af are the dipole-dipole interaction strength as in Eq.(19). 
Therefore, nuclear spins now evolve independently of each other and couple with the NV center individually, see 
Eq.(27). 

Polarization dynamics of the nuclear spins. — The whole polarization process consists of repetitive cycles. In each 
cycle, we initialize NV center in the \ms = 0) = |^) state with a green laser (537 nm), and prepare it in the state 

1^^) ~ ^ 1^)) ^^^^ ^ f ™crowave pulse. After one polarization cycle, the nuclear spin state evolve as 

Pf{nr + r) = trnv[exp {-itHp)p^{nr) exp {itHp)] (28) 

where r is the time of each cycle. It can be seen from the effective Hamiltonian in Eq.(27) that, when the effective 
spin-^ is on resonance with nuclear spins, polarization can be efficiently transferred from the NV center to the nuclear 
spins. The appropriate resonant condition, called Hartman-Hahn condition [S7], turns out to be 

uJnv = ^nv - (7ivB - Ap) =ujf = ^A^ + (l]p/2)2. (29) 

The dynamical nuclear polarization process will finally polarize nuclear spins towards the product state |i) (8) • • • (8) 
We have performed exact numerical simulations for a 3x3 lattice with a Chebyshev expansion to calculate 
U{r) = exp {—itHp) where Hp is the Hamiltonian in Eq.(20-22) with the rotating wave approximation, by assuming 
that the magnetic field is much stronger than the nuclear spin coupling, which is easily satisfied in experiments. 
The result is shown in Fig.3(a-b) of the main text, which shows that nuclear spins are indeed isolated and our 
polarization scheme works very well for the present system. Due to the limited computational power, to estimate 
the polarization efficiencies in large systems, we describe the polarization dynamics with the following master equation 



-Pf = - - -(pfA+A_ + A+A_pf-2A_pfA+) (30) 

where = Zli^'l^f ^± = ^iOt^f- The above master equation can be derived from Eq.(28) by assuming that 
the polarization cycle time r is sufficiently short (namely satisfies T^^^{g^Y = e <C 1) [S8]. To numerically solve 
the above master equation, one needs to make some approximations on the coherence between nuclear spins. If we 
neglect coherences between nuclear spins (known as the spin temperature approximation), then the polarization rate 

is Pi = {g^Yr = ^{gj^Y / ^^]^{g^Y ^ ^{eg^). In this case, the polarization time scales linearly with the total 
number of nuclear spins. To take into account the coherence between nuclear spins, one can approximate the higher 
correlation terms by a Wick- type factorization [S8] as follows 

i{s+s|s-) = i(s+s-) - {s+s-)(s+s-) + (s+s-){s+s-). (31) 

Comparing these approximations with each other and with the exact simulations on small systems can help us to 
obtain reliable estimations for the polarization efficiencies in large systems. In Fig3.(b) of the main text, we compare 
the exact numerical calculation with the approximation of spin temperature for a 3 x 3 nuclear spin lattice. It can 
been seen that for the chosen r, they show good agreement. We also compare the polarization dynamics for a 10x10 
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t/N (ms) 



FIG. 4: The average probability of nuclear spin down state as a function of the polarization time r. Comparison of the 
spin temperature (red, diamond), Wick-type (green, circle) for a 10 x 10 nuclear spin lattice with the polarization cycle time 
T = 3/is. The spin-temperature limited polarization efficiency can be approached by introducing magnetic noise to remove 
nuclear spin coherence every 20ms, the result is calculated with the Wick-type approximation (purple, square) with r = 

The magnetic field direction is m = (\/ 1? 0? 



nuclear spin lattice under the spin temperature and the Wick-type approximation as shown in Fig. 4. Thus it is 
possible to estimate the polarization rate based on the spin temperature approximation, i.e. the polarization rate is 
Pi ^ j^i^di')- From these results we estimate that for the example of a fluorine nuclear spin array with a distance 
of 5nm from the NV center, the polarization time scale for a total number of N nuclear spins is {pi)~^ = N • lO(ms) 
with e = 0.1. Further improvements of the polarization efficiency can be expected by optimizing the cycle time. 

NV center as a quantum probe for measurement of the nuclear spin state. — Following initialization and evolution 
of our nuclear spin quantum simulator, we must be able to determine its properties via appropriate measurements. 
The measurement of individual nuclear spins is challenging due to their small magnetic moment. In the main text, we 
propose that an NV center can serve as a measurement interface for nuclear spin states. We adopt the same strategy 
for isolating nuclear spins from each other during measurements as in the polarization process described in the main 
text. Thus, we need to apply an RF-pulse to map the nuclear spin states from the original spin basis |t)} to the 
basis for measurement The dynamics is again described by the Hamiltonian in Eq.(27) 

= + -f E ^1 + E s!^' + E 9^ Hf^ + ^•^) • (32) 

i i i 

Before measurement, we prepare the NV center either in the state It^^) or ||^). The microwave driving field on the 
NV center is then tuned to match the Hartmann-Hahn condition, see Eq.(29). After time r, we perform measurement 
on the NV center which is denoted as an observable A, which can be written as follows 

Air) = AiO) + i^[Hp,A] - ^[Hp, [Hp, A]] + ... (33) 

In particular, we can measure the population of states |+) = and |— ) = \lx) of the NV center, and the corre- 
sponding observables are A = |+)(+| and A = \—){—\ respectively. One can obtain 

P-= r'Y.Y.^9i9f){n~^-) = r'Y.^9i?P} E iot 9f)(^^ + K^) ^ (34) 

P+= ^'T.T.(9ht)iK~sp=^'T.(9h'P^ + r' E (ghtMsJ+s-sl) (35) 

* 3 * ih3)^7^3 

where Pj and P^^ are the spin up and down population of individual nuclear spins, j) i^j over all pairs 

of nuclear spins. Therefore, we obtain 

Ap = P+ - = Y.^gmPl - P}) = Y^^difiK)- (36) 
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The observable Ap thus provides information about the average magnetization of nuclear spins. We also obtain 

AxY = Pl + P+ = E T.(9ht){^nj + = E Y.(9^9t){^r^- + sfsj) (37) 

i j i j 

which provides information about {xx and yy) correlation functions of nuclear spins. To estimate the other correlation 
functions, we introduce the Hadamard operation Oh^ the phase transformation Oj and Op as follows 

The spin operators under the transformations Oh and 0/ are listed as follows: 

Oj^s^O// = s^ OljsyOH = -s^ Oj^s^O/f = s^, (39) 
0}s^0/ = -s^ 0}s^0/ = -s^, 0\s'Oi = s^ (40) 
Oj^s^Op = s^ Oj^s^Op = -s^, Ol^s'Op = -sy. (41) 

Therefore, we apply an RF-pulse corresponding to the transformation O = Oh^Op on the nuclear spin state |^) 
0|^) such that the measurement of the observable, as described in Eq.(37), provides information about the correlation 
functions (yy + zz) or {xx + zz) of the fluorine nuclear spins. 

Ayzli*) = Axy|o„|*) = 2r2EE ((^L ® O^) (§f §1 + {Oh^E^Oh)) (42) 

* j 

= 2r2 E E ialaf) (sf + sf sj) , (43) 
Axzli*) = A^yloH*) = 2r2EE ((4 ® 4) (Sfs,^ + §i §i) (Op ® Op)) (44) 

* j 

= 2r^EE(5.V)(§r§.^ + (45) 

« j 

Together with Axy in Eq.(37), it is thus possible to estimate all the correlations of XX, YY and ZZ of nuclear spin 
states. 

To measure the observable such as nuclear structure factors of nuclear spin state, we can apply a gradient field on 
the nuclear spins. Accordingly, the nuclear spin at the position experiences a field bi = Vi ■ (hx^hy)^ and gains a 
position-dependent phase (/;^ = • q where q = {qx^Qy) = {jNtp) {bx^ by). We then perform the same measurement as 
in the above discussions and obtain 

E (aht) cos m - ?i) • q] i^r^j + sfsj) , (46) 
E (aht) cos m - ?.) • q] (§r§I + , (47) 
E (aht) cos m - ?i) • q] (§f §i + • (48) 

* j 

This makes it possible to extract the information on the structure factor defined as follows 

5q {q.,qy) ^ ^ ^1 E^'^^'^fel') «^ (aL + Ail) + (a^^ + A^l) - (a|^ + A-^l) . (49) 

The accuracy of the above estimation depends on the difference of the couplings between the NV center and 
individual nuclear spins. For the relevant nuclei (e.g. with non-zero correlation functions), their mutual distances 
are usually much smaller than their distances from the NV center. Therefore, their individual couplings with NV 
center (i.e. g:^ and gj- in Eqs. (46-48)) will be similar. For translational invariant systems, the observable estimated 
from NV measurement is expected to give us information on the average properties of the system state. As the 







= 2t' 


A^z- 


^ A"^ 


= 2r' 


A|z- 


h A"^ 
^ ^xz 


= 2r' 
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observable is estimated by the measured quantity P^/r'^ with /i, = ±, the measurement accuracy wih depend on the 
choice of measurement time r. In practice, one can perform measurements for different values of r and estimate the 
observables with improved accuracy. In the situation of quantum phase transitions, which are usually accompanied 
by sudden or non-analytic changes of observables, the measurement via NV centers is expected to be able to identify 
accurately different phases, e.g. see Fig. 5 in the main text. 

Tuning parameters for quantum magnetism. — The nuclear spin coupling strength can be controlled by the magnetic 
field direction as 

ai-^i) = ^ - 3 (r.. ■ zf] = ^ (1 - 3cos^ e,,) (50) 



47rr3. 



where v^i = rkir%i is the vector that connects the nuclear spin k and I, and Oki is the relative angle between r^; and 
the magnetic field direction z. In the main text, we consider two examples of magnetic fields to tune the relative ratio 
between the nearest neighbor nuclear spin interactions along the three directions, i.e. Jai, Jo2j J&s with ai = (1,0,0) 

a2 = (i^,0),a3 = (-^,^,0). 




2/3 1 2/3 1 0.2 0.4 0.6 0.8 1.1 0.2 0.4 0.6 0.8 1.0 

cos^ cos^ cos^ cos6' 



FIG. 5: The nearest neighbor interaction strength Ja^ in the unit of ga — 6.8kHz as a function of the direction of magnetic 
field rfi — (0, cos^, sin^) (a) and rh — (cos^, 0, sin^) (b): Ja^ = Ji (red), = Jd^ = J2 (green). The ratio = = as 
a function of cos^ for the magnetic field direction m = (0, cos^, sin^) (c) and m = (cos^, 0, sin^) (d). 

In the first example, where the magnetic field direction is m = (0, cos6>, sin6>), the nearest-neighbor couplings 
are Jq,^ = Ji = g^ (AF), Ja^ = Ja^ = J2 = 9a{^ — |cos^6>), together with the ratio are plot in Fig.5(a-b). For 
cos^ > |, we obtain J2 < and thereby are ferromagnetic interaction. For the particular value cos^ = |, J2 = 0, 
and the system consists of ID (AF) chains interacting through weaker non- nearest-neighbor couplings. In the other 
example, where the magnetic field direction is m = (cos 6>, 0, sin6>), the nearest-neighbour interaction strengths are 
Ja^ = Jl = g{a){l — 3cos^ 6>), Ja^ = Jds = J2 = — | cos^ 6>), as represented in Fig.5(c-d). As the magnetic field 

direction changes, J2 is always positive and thus is antiferromagnetic (AF) , while Ji is negative (ferromagnetic) for 
cos l9 > 1/3. 

Comparison of quantum magnetic phases with the short-range model. — In the main text, see Fig. 5, we present 
two examples of quantum magnetic phase transitions that are controlled by the orientation of the external magnetic 
field. Here, we have taken into account long-range interactions up to the cut-off distance of >/T3a o:^ 3.6a, where 
a is the lattice constant . In the first example with the magnetic field direction m = (0, cos^, sin^), we find that 
the competition between nearest-neighbor interactions Ji and J2, together with the strongest next-nearest neighbor 
interaction, leads to a transition towards the ferromagnetic phase as the value of cos^ increases, see Fig. 5(a) in 
the main text. We remark that the non-nearest neighbor interaction in fact plays an essential role in such a 
phase transition. For comparison, we calculate the spin structure factors for the short-range model with only 
nearest-neighbor interactions. It can be seen from Fig. 6(a) that the ferromagnetic phase does not appear in the 
corresponding parameter regime. 

Let us move to the second example, where the magnetic field direction is m = (cos 6>, 0, sin6>). By only considering 
nearest-neighbor interactions (which are all negative values), the system does not display frustration, and we obtain 
the ferromagnetic order in the ai direction, and the anti- ferromagnetic order in the a2 direction (F-AF phase), 
see the corresponding spin structure factors for this short-range model in Fig. 6(b). Moreover, as the value of cos^ 
becomes smaller, J2 decreases while the non-nearest-neighbor interactions now becomes comparable with J2, see 
Fig. 5(b) in the main text. Then, the triangle consists of J2 and J31, J32, J33 (see Fig. 5(a) in the main text) becomes 
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FIG. 6: Quantum magnetic phase transitions of fluorine quantum simulation on a triangular lattice for short-range models, i.e. 
only nearest-neighbor interactions are taken into account, (a) The spin structure factors changes with the magnetic field in the 
direction rh — (0, cos^, sin^) as compared with Fig. 5(a) in the main text, (b) The spin structure factors changes the magnetic 
field in the direction m = (cos^, 0, sin^) as compared with Fig. 5(b) in the main text. The spin structure factors are calculated 
using the Lanczos algorithm for a 6x4 lattice under periodic boundary condition. 



frustrated, and the competition between the nearest-neighbor interaction and long-range interactions leads to a new 
magnetically order phase. The spins are ferromagnetic in the ai direction, while anti-ferromagnetic in the sublattice 
every second line in the a2 direction, see Fig. 5(b) in the main text. 

Quantum magnetic phase transitions with fluorine defects. — We have demonstrated different magnetic phases in 
Fig. 5 of the main text. Here, we demonstrate that key signatures for these phase transitions persist even when the 
lattice has defects, e.g. the dangling bond of carbon on diamond surface is terminated by spin-less nuclei rather than 
fluorine. We calculate the spin structure factors under the same magnetic field condition as in the main text, on a 6x4 
lattice with more than 20% sites at random positions being absent. The results, averaging over 30 random realization 
of the triangular lattice with defects, are shown in Fig. 7, which provide evidence that the transition between different 
quantum magnetism phases are still observable in the presence of significant disorder. We remark that the exact 
diagonalization is only performed on a relatively small lattice due to the limited computational overhead. Further nu- 
merical simulations will be necessary to determine the exact quantum phase diagram for various levels of defects rates. 




-1 -0.5 -8 -7 -6 -5 -4 -3 -2-10 1 

J2/J1 J2/J1 



FIG. 7: Quantum magnetic phase transitions of quantum simulation on a triangular lattice with defects, (a) Tune the phase from 
ID anti-ferromagnetic chains to ferromagnetic with the magnetic field in the direction m — (0, cos^, sin^) on a 6x4 lattice with 
25% site depletion, (b) The phase transition from ferromagnetic-antiferromagnetic (F-AF) order to ferromagnetic- (alternative) 
antiferromagnetic (F-NAF) controlled by the magnetic field in the direction m — (cos^, 0, sin^) on a 6x4 lattice with 20% site 
depletion. The spin structure factors are calculated using the Lanczos algorithm under periodic boundary condition. 

Hard-core boson superfluid and supersolid. — The nuclear spin Hamiltonian of our system can be mapped to the 
hard-core boson model as 

5 = ^ (VijUiTij - tij (^alaj + ttitt]^^ + /i ^ ni (51) 
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FIG. 8: Enhancement of superfluidity by long-range interaction, (a) Comparison of superfluidity between our model with 
long-range interaction and the short-range model with the system size 24 x 24. (b) Robustness of superfluidity with long- 
range interaction under the depletion of lattice with the system size 20 x 20: the site depletion probability is 25% (circle) and 
50% (square). The geometry frustration (as characterized by the ratio V2/V1) is controlled by the magnetic field direction as 
m = cos^as + sin^Z with Z = (0, 0, 1). The temperature used in QMC simulation is T/ga = 0.1. 



by the Holstein-Primakoff transformation [S9] = — ^ (n^ = 0,1), = sf -\- iSy = ajv^l — rii and 
s~ = sf — iSy = ^/Y^^riiai. Here, the chemical potential is /i = {jnB) — ^jSi'^ij)^ the repulsive interaction is 
Vij = gi^ij)^ and the hopping rate is tij = ^g{Yij). Depending on the ratio between the tunneling, the repulsive 
interaction y = ^ and the lattice geometric frustration, the system may demonstrate interesting phases such as the 
solid (S), superfluid (SF) or supersolid (SS) phases. The superfluid is a phase which possesses long-range off-diagonal 
order; while the supersolid exhibits both long-range off-diagonal and diagonal order, showing thus the features of 
a superfluid and a solid simultaneously. Our model has a large flexibility in tuning the geometric frustration and 
the ratio between the repulsive interaction and the hopping, both of which are long-range. It is thus possible to 
investigate rich phases of hard-core boson models. On the other hand, the model with frustration and long-range 
interactions poses great challenges to the classical numerical simulations. Quantum Monte Carlo (QMC) simulation 
which might be efficient for 2D systems are likely to become inefficient due to the sign problem. Here, we present sim- 
ple examples which nevertheless can provide interesting insights into the properties of superfluid and supersolid phases. 




FIG. 9: The superfluid density ps and the normalized structure factor see the illustration in the inset of a) as 

a function of chemical potential /i — /io, where /io is the chemical poterntial corresponding to half-fllling, for a lattice of size 
16 X 16 (a) and 20 x 20 (a) with the parameter t/V = 0.2. We also plot the curves for the results of a lattice of size 12 x 12 

for visibility. The magnetic fleld direction is ?fi = cos^as + sin^Z with cos^ = and Z — (0, 0, 1). The temperature used in 

QMC simulation is T/ga = 0.1. 

With the magnetic field along the direction rh = cosOas + sin^Z, where Z = (0,0,1), the nearest-neighbor 
interactions become Va^ = = Vi = ga{l — |cos^6>), and Va^ = V2 = ^a(l — 3cos^6>). By changing the value of 
the magnetic direction angle 6>, we can gradually tune the geometric frustration, as quantified by the ratio V2/V1. 
The geometric frustration is important for the emergence of the superfiuid, as can be seen from the result for the 
short-range model shown in Fig. 8(a). It can also be seen that taking long-range interactions into account will signifi- 
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FIG. 10: Quantum phase transition of superfluid (SF) and supersolid (SS) with 10% lattice depletion, (a) The average filling 
(n) as a function of chemical potential /i — /xo, where /jq is the chemical potential corresponding to half- filling. We also plot the 
curve for the same model without lattice depletion for visibility (square). The kink is indicated by the arrow, (b) The superfluid 
density ps and the normalized structure factor Sg*(7r, tt/v^)- The system size is 12 x 12 (blue, circle), 16 x 16 (green, diamond) 

and 20 x 20 (red, triangle). The magnetic field direction is m = cos^as + sin ^(0,0, 1) with cos^ = ^J^. The parameters are 
the same as Fig. 6 in the main text: t/V = 0.2, and the temperature used in our simulation is T/ga = 0.1. 



cantly enhance the superfluidity, which is robust under a relatively high level of lattice depletion (defects), see Fig. 8(b). 

If we choose cosO = for the magnetic field direction m = cosOas + sin^Z, the lattice is equivalent to a 2D 

rectangular lattice. The model incorporates long-range repulsive and hopping interactions. For the ratio t/V = 0.2, 
we have calculated the superfluid density ps and the normalized structure factor Sg*(7r, tt/a/S) for a lattice of size 
12 X 12, which evidences the solid (S), superfluid (SF) and supersolid (SS) phases, see Fig. 6 in the main text. In Fig. 9, 
we show the results for a lattice of size 16 x 16 and 20 x 20, and find that both ps and Sg'(7r, 7r/\/3) are finite and 
size independent in the regime fi — fio e [0.6,0.8]. This provides evidence supporting the existence of the supersolid 
phase. We have also performed simulations including the lattice site depletion. The results are shown in Fig. 10. As 
compared with Fig. 6 in the main text, the lattice depletion would downgrade the superfluidity. Nevertheless, the 
supersolid phase may still exist with a certain level of lattice depletion (10% in our calculation). The full study of 
quantum phases with lattice depletion for the present quantum simulation will be very interesting, which is beyond 
the scope of the present work and will be addressed in detail in a future publication. 



Quantum quench and spin squeezing. — To illustrate our ideas, we consider the dynamics of a fluorine quantum 
simulator on a triangular lattice which undergoes a quantum quench. After the polarization, the system is prepared 
in the product state |'^(0)) = (g) ••• (g) which is the eigenstate of the Hamiltonian with a large resonant 
RF-field Hsq = j) di^ij) [^i^j ~ ^ (sf s| + ^j)] + ^sq Xli • Then we suddenly decrease ftsq to a value that is 
comparable with g{Tij) (in our numerical calculation, we use the value of lOkHz). After such a quantum quench, the 
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FIG. 11: Squeezing parameter (a) and (b) as a function of the evolution time t for the fluorine simulation on a 4x4 
triangular lattice. We compare the estimation with NV center measurement (purple square) with the exact values (green 
diamond) . 
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state l'^(O)) is not anymore the eigenstate of the new Hamiltonian and the nuclear spins start to entangle with each 
other. We remark that during the dynamics, the NV center is polarized to the state = to avoid its effect on 
the nuclear spins. We characterise the system dynamics by the quantity of spin squeezing, which shows a reduced 
variance of the collective angular momentum even when the mean angular momentum in an orthogonal direction is 
large. The squeezing parameter is = 2 (AJc^)^ /|(Ja;)|, where the collective angular momentum is = ^j^I^ 
and (AJ«)' = (J2) - (J^)2 with a = x^y^ z. Spin states with < 1 are referred to as spin squeezed states and 
necessarily entangled [SIO]. We perform exact numerical simulation for a 4x4 triangular lattice and calculate the 
evolution of squeezing parameters, see Fig. 11. This shows that the angular momentum variance in the y direction 
is reduced. To measure the squeezing parameters with the NV center, we apply the Hadamard operation Oh and 
then measure the nuclear spin magnetization which gives us the estimation of (Jx)- The value of (J^) can be 
approximated by {A^is + Aq,^ — A/3^)/4. For the simplicity of numerical calculation, we choose the measurement time 
r = min{P^|T- = 10%,42/i5}, and the observable is estimated by the measured quantity P^/t'^ with /i, = ±. 
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